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Abstract 

A numerical investigation of grain-boundary (GB) grooving by means of the Level 
Set (LS) method is carried out. GB grooving is emerging as a key element of electromi- 
gration drift in polycrystalline microelectronic interconnects, as evidenced by a number 
of recent studies. The purpose of the present study is to provide an efficient numer- 
ical simulation, allowing a parametric study of the effect of key physical parameters 
(GB and surface diffusivities, grain size, current density, etc) on the electromigra- 
tion drift velocity as well as on the morphology of the affected regions. An idealized 
polycrystalline interconnect which consists of grains separated by parallel GBs aligned 
normal to the average orientation of interconnects surface is considered. Surface and 
grain-boundary diffusion are the only diffusion mechanisms assumed. The diffusion is 
driven by surface curvature gradients and by an externally applied electric field. The 
corresponding mathematical system is an initial boundary value problem for a two- 
dimensional Hamilton-Jacobi type equation. To solve for the electrostatic problem at 
a given time step, a full model based on the solution of Laplace's equation for the 
electric potential is employed. The resulting set of linear algebraic equations (from 
the finite difference discretization of the equation) is solved with an effective multigrid 
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iterative procedure. The details of transient slit and ridge formation processes are pre- 
sented and compared with theoretical predictions on steady-state grooving pi 0, 0]. 

Keywords: Level Set method, modeling, electromigration, grain boundary grooving, 
drift. 

1 Introduction 

This paper is a continuation of our work on numerical modeling of the formation and 
propagation of groove-like defects at GBs in thin film polycrystalline interconnects used 
in microelectronics (ME). 

In modern ME industry, the reliability of ME integrated circuits has become no 
less important than their performance. Some of the most vulnerable elements of ME 
circuits, susceptible to several types of failures, are the interconnects. These are thin 
film metallic conductors which connect the active elements. 

The defects (due to the small cross-section, high current density, mechanical stresses 
and presence of GBs acting as fast diffusion pathways) lead to a loss (in relatively 
short times) of electrical and mechanical integrity, i.e. to line opens or shorts. For 
example, in the presence of a large GB flux (J g b « 10 -4 iim 2 /s) a groove can extend 
several micrometers in a few hours [p|. Thus, GB grooving is one of the main failure 
mechanisms in advanced integrated circuits. 

In the absence of an external potential field and mechanical stresses, the GB atomic 
flux Jgj, = 0, and the corresponding groove profile evolves via surface diffusion under 
well-known conditions of scale and temperature (the so-called Mullins problem Q ) . In 
this case, mass transport by surface diffusion is driven only by the surface Laplacian 
of curvature. Essentially, matter flows from low-curvature regions to high-curvature 
regions. 

In ||, we presented and discussed the numerical approach (e.g. the Level Set 
(LS) method) used to model GB grooving phenomena. We also tested the LS method 
on two simple, already solved, grooving problems: Mullins problem, and that of GB 
grooving by surface diffusion in a periodic array of stationary GBs Q. In both cases, 
the results obtained by means of the LS method are in good agreement with the 
theoretical predictions. In this paper, we consider the second geometry only, as being 
more realistic (see Fig. [l|). Due to axial symmetry at x = 0, x = L (where L is the 
grain size), we do not attempt to calculate groove branches at x < 0,x > L. 

Electric fields/currents in metallic conductors provide an additional driving force 
for surface/GB mass fluxes 0. In the presence of an electric field, collisions between 
the conduction electrons and the metal ions lead to drift of the ions. This process is 
known as electromigration (EM). 

GB grooving with a GB flux in real thin film interconnects is a complex problem. An 
adequate numerical modeling technique should be capable to manage such issues as GB 
grooving with an arbitrary EM flux, and various ratios of GB to surface diffusivities; 
the latter was predicted to critically affect groove kinetics and shape, and thus account 
for various EM failure regimes (see jj], ||, [3| and the references therein). In cited 
works, analytical and semi-analytical approaches for analysing steady-state grooving 
regimes were employed. However, to our knowledge, no effort has been made to directly 
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numerically simulate the transient stage during GB grooving. This is the ultimate goal 
of this paper. 

We do not consider mechanical stresses in GBs which, as a matter of fact, are invari- 
ably induced by the field |7| (the approximations under which it is reasonable to neglect 
the stress are discussed in [||, [|). Also, under typical operational conditions of ME 
interconnects, lattice transport may be neglected compared to surface/GB transport 

§!• 

Our paper proceeds as follows. In Section 2, we give details of the physical formula- 
tion. In Section 3, we discuss some improvements in the numerical algorithm and also 
the aspects which are due to incorporation of the electric field/GB flux in the model 
(for other algorithmic details, the interested reader should refer to ||). Our numerical 
results and discussion are presented in Section 4. 



2 Physical model 

2.1 Driving forces for the diffusion 

In the absence of an electric current, the surface diffusion is driven by a variation 
in chemical potential, fj, s , which causes atoms to migrate from high potential to low 
potential regions. It was shown M that 

H S = K sls n (2.1) 

where K s is the surface curvature, "y s is the surface tension, and O is the atomic volume. 
Gradients of chemical potential are therefore associated with gradients of curvature. 

Let r be the tangential direction to the surface profile in 2D. Let x, y be Cartesian 
coordinates along horizontal and vertical boundaries of the computational box (Fig. 
|l]). If n = (n x ,n y ) is the unit vector normal to the surface, then the following relations 
hold: 

8K S _^ 8K S 8K S 

— — = \/K s ■ t = -^—n y — 

or ox ay 

The corresponding surface flux (volume crossing unit length per unit time) is then 
given by 

J ° -~ltr-fr-- BK T (2 ' 3) 

where the superscript indicates that the flux is due to the curvature gradient, 

B= D i Sgl (24) 

is known as Mullins constant, and D s , S s , k, and T denote surface diffusion coefficient, 
thickness of the surface diffusion layer, Bolzmann's constant and absolute temperature, 
respectively. Note that J^ Ks is proportional to the first directional derivative of the 
curvature. 

If an electric field is present, the flux J s of matter at the curved surface of the 
conductor is driven simultaneously by curvature gradients, and by the component E of 
the local electric field along the surface. In what follows, we distinguish between two 
models which handle the electric field: 



(n y , -n x ), = VK S ■ t = -^-n y ^-n x = K S T . (2.2) 
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1. Piecewise-constant electric field 

Let C and O denote the conductor (interconnect) material domain and the outer 
(surrounding) material domain above the surface profile, respectively (see Fig. 
[l]). In this model, the vector of the electric field intensity is parallel to the GBs, 
Eo = (0, Eq v ), and Eq v is a step function, 

, . _ j E l n = const. if grid point (x^yj) eC . . 

We assume that the surrounding material is less conductive than the interconnect 
material, and therefore I-Eq"'! < I -So™ I- ^ n our numerical experiments we chose 
the ratio I-Eq^I/I-Eq"! = 0-1- in the finite difference approach, the discontinuous 
distribution of electric field intensity is smoothed out across the surface profile 
(see the details in Section 3). The component of the local electric field along the 
surface, E, is then approximated by the projection of Eo on the surface, E = Eo-t. 
This approximates the true value given by solving Laplace's equation for the 
potential, subject to the boundary conditions of constant fields of magnitudes 
Eq 1 and Eq u1 in the conductor and surrounding material domains, respectively. 
The corresponding electrically induced surface flux of matter is given by 

Jf = -£^±E = -B e E (2.6) 
where the superscript indicates that the flux is due to the electric field, and 

Be = ^ (2.7) 

where Z s = z*e is the effective charge of the ions undergoing electromigration 
in the surface layer and e is the unit electronic charge; the sign of z s is usually 
positive (i.e., matter flux in the direction of the electron flow). 

2. Solution of Laplace's equation for the potential 

Assume that (at a given time step of overall marching algorithm) U(x,y) is 
the electric potential within the (rectangular) computational box. U~{U + ) and 
U + (U~) are its values on the upper and lower boundaries of the box, and U n 
is the normal derivative on the boundary. U~ and U + are assumed to be time- 
independent and uniform along the boundaries; U + — U~ is the external voltage 
applied to the interconnect. The distribution U(x,y) is governed by a static 
elliptic PDE 

d /, dU\ d /, dU 



dx \^ dx ) dy \ dy ) ^' ^ ^ 

with boundary conditions U n = U x = on the vertical boundaries of the box 
(which in our case coincide with GBs). Equation ( |2.8] ) is derived from the well- 
posed three dimensional potential problem for the two-layer interconnect. The 
assumptions and complete derivation for the case of small aspect ratio are pre- 
sented in [|10f| . We also give some details in the Appendix. In eq. (p.8[ ), k = k(x, y) 
is the specific electrical conductivity (at a given time step) of the material which 
fills the computational box. To solve (|2.8| ), a finite difference scheme was devel- 
oped and analysed in (1^]. The distribution of the specific conductivity in the 
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physical system under consideration is discontinuous: the conductivity inside the 
conductor material (domain C, Fig. |l]) differs by a finite value from that of the 
surrounding material (domain O). We assume 

k _ j hn = const. > if grid point (xi, yj) € C , , 

1 k out = const. > if grid point (xi,yj) € O, 

i.e. k = k(y) is a step function. In our numerical experiments we chose the ratio 
kout /kin = 0.1. Since the surface of the conductor evolves in time and space, 
then to find the time-dependent solution U(x, y, t) we need to solve the static 
equation (|2.8|) every time step with k given by fl2.9|) . In order to be able to 
compute accurately the electric field intensity (which is the derivative of U), the 
discontinuous distribution of the specific conductivity is smoothed out across the 
surface profile. The finite difference discretization of ( [2.8[ ) in the computational 
domain leads to a set of linear algebraic equations with a sparse banded matrix. 
This set is solved with an effective multigrid iterative procedure |l(| . The solution 
of the previous time step is used as an initial approximation for the current step 
which allows fast convergence. 

After the potential is established everywhere in the computational domain, the 
corresponding electrically induced surface flux jf is given by (pH|), where 

E = -t-VU (2.10) 



To summarize the above discussion, the total flux of matter along the surface is 

J s = J™ + Jf. (2.11) 



Physically, equation ( 2.1l| ) says that atoms will diffuse in the direction of the electron 



flow if the field dominates, but toward the position with the large curvature if the 
surface energy dominates. This competition between the electric field and the surface 
energy is essential for the groove dynamics. 

The electric field results also in the diffusion of matter along GBs. The diffusion 
flux along the GB, J g h, is given by 

J 9b = -M^E (2.12) 

where Dgb, 5 g b, Z g b = z* b e > are the GB diffusion coefficient, thickness and effective 
ionic charge, respectively, and E is the component of the electric field along the GB. 



2.2 Boundary conditions at groove roots 

The evolution of the surface is constrained by two conditions imposed locally at groove 
roots a and b (Fig. Q): 

1. Equilibrium angles 

The boundary condition is dictated by the local equilibrium between the surface 
tension, j s , and the GB tension, y g j,. In the symmetric case of a GB (x = 0) 
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normal to an original (y = const.) flat surface, the angle of inclination of the 
right branch of the surface at the groove root with respect to the x axis is M\ 

6*o = sin~ 1 ( r y g i ) /2 r y s ) = const. (2-13) 

The rapid establishment of the equilibrium angles between the GBs and the sur- 
face by atomic migration in the vicinity of the intersections develops some curva- 
ture gradients at the adjacent surface, and thus induces surface diffusion fluxes, 
J^ K , along the groove walls. The directions of the fluxes depend on the sign of 
the respective surface curvature gradients at the groove groots. 

2. Continuity of electrically induced fluxes 

The boundary conditions read 

J gb = 2jf , (2.14) 
since both branches of the groove act as sinks or sources of matter. 



u~ (u + ) 

GB, outer material domain (O) 




interconnect material domain (C) 



U + (U~) x L 

Figure 1: Sketch of GB grooving in a periodic array of stationary GBs. The grain size is L, 
and groove root points are marked as a and b. 



3 The numerical procedure 

The Level Set method is used to "capture" the evolution of the conductor surface. The 
method was introduced by Osher and Sethian [11] and was further developed during 
the last several years. The method enables to capture drastic changes in the shape of 
the curves (surfaces or interfaces) and even topology changes. 
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The basic idea of the method consists of embedding the curve y(x, t) into a higher 
dimensional space. As a matter of fact, we consider the evolution of a two-dimensional 
field (j)(x, y, t) such that its zero level set, cj>(x, y, t) = 0, coincides with the curve of 
interest, y(x,t), at any time t. The level set function c/>(x,y,t) can be interpreted as a 
signed distance from the curve y(x,t), which moves in the direction normal to itself. 

The evolution of 4>(x,y,t) is described by an Hamilton-Jacobi type equation. A 
remarkable trait of the method is that the function <j>(x, y, t) remains smooth, while 
the level surface = may change topology, break, merge, and form sharp corners 
as (p evolves. Thus, it is possible to perform numerical simulations on a discrete grid 
in the spatial domain, and substitute a finite difference approximations for the spatial 
and temporal derivatives in time and space. Another nice feature of the method is 
that the explicit location of the interface needs not to be known in the computational 
process; all the necessary information is extracted from the level set function. 

The evolution equation has the form 

4 + F|V0|=O given cf)(x,t = 0). (3.1) 

The normal velocity, F, is considered to be a function of spatial derivatives of 4>(x, y, t). 
In many applications F is a function of the curvature, K s , and its spatial derivatives. 
The curvature K s may be computed via the level set function as follows: 



K s = V • n, n 



|V0| 



(3.2) 



Here n is a "normal vector", and it coincides with the (previously introduced) unit 
normal to the surface, y(x,t), on the zero level set = 0. Formulas ( |3.2[) can be 
combined as follows 

v „ V0 <Pxx4> 2 y ~ 1<\>x4>y4>xy + <f>yy4>l fo x 

K -^-wr — — • (3 ' 3) 

and the sign of K s is chosen such that a sphere has a positive mean curvature equal to 
its radius. In the case of surface diffusion in 2D, 

where J s is given by fl2.11| ). 

The difficulties in the numerical solution of ( |3.1[ ) in our case are due to the fact 
that, as could be noted from ( |2.3[ ), ( |3.3P , fl3.4| ), the first term in F contains space 
derivatives of order 4 of the level set function. Therefore, the evolution equation (|3.l| ) 
is highly sensitive to errors. Besides, this fourth derivative term leads to schemes with 
very small time steps. 

In ||, we presented the computational algorithm which solves the problem of GB 
grooving by surface diffusion in the absence of electromigration. This could be viewed 
as the limiting case of the problem which is under consideration in this paper, corre- 
sponding to the situation where electrically induced surface and GBs fluxes jf and 



Jgi, vanish. The normal velocity function (3^) in the latter case contains only the first 



term. The basic features of the algorithm are: 
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• the use of a uniform static grid in both space directions 

• the use of a standard second order-accurate finite difference scheme in space 

• the approximation of spatial derivatives (in normal direction) on the boundaries 
of the computational box by second order one-sided differences 

• time marching is done by a second-order Total Variation Diminishing (TVD) 
Runge-Kutta procedure jL2L O 



the use of second-order Essentially Non-Oscillatory (ENO) scheme [14] to approx- 
imate the gradient function in ( p.l| ) 

the use of "reinitialization" |l5| every time step to keep the level set function <f> 
a signed distance function 



The solution of ( |3.1|) (in Mullins case of an infinite bicrystal with a single GB) subject 
to the conditions of a constant angle of surface inclination and zero surface flux J^ K 
at the groove root a, is then a self-similar surface profile, whose linear dimensions are 
proportional to (Bt) l / A , B given by ( |2.4| ) (Fig. |2|). If the dimensions of the crystal 
are finite, grooves develop at each GB; grooving stops when, at sufficiently long times, 
identical circular arcs develop connecting adjacent GBs (Fig. ^). The parameters 
chosen for the runs are typical for copper interconnects at temperatures relevant to 
experiments (T = 600 K) [|, ge]]: Q = 1.18 x 10~ 29 m 3 , D s = 3.3 x 10" 14 m 2 /s, j s = 
1.7 J/m 2 , 5 S = 3.5 x 10~ 10 m, kT = 8.28 x 10~ 21 J. It is worth noting that our numeri- 
cal treatment of GB grooving is not constrained by the assumption of small equilibrium 
angles ( "small slope approximation" ) , in contrast to the analytical approaches of the 
pioneer works [||, f|. 
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0.36 




0.3 1 ' ' ' ' 

0.5 1 1.5 2 

x, |xm 



Figure 2: Groove development and propagation along the GB of an infinite bicrystal (Mullins 
problem). Grain size L = 2 /an. For copper interconnects at T = 600 K (327 °C), B = 
9.2 x 10~ 33 m 4 /s. The angle at the groove root 6 = tt/22 (tan# « 0.144). 
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Figure 3: GB grooving in a periodic array of stationary GBs. Parameters are the same as 
in Fig. H 



In special attention was given to the treatment of constant-angle and zero-flux 
conditions at the groove roots within the framework of the Level Set method. Two 
methods were developed, the first based on interface reconstruction every time step, 
with subsequent correction of the angles followed with reinitialization, and the second 
based on the extension of the 0-field beyond the GBs using the expansion in Taylor 
series up to second order. Both methods were successfully used in calculations, but we 
observed that sometimes both procedures resulted in a loss of accuracy. In this paper, 
we propose a new, robust and highly accurate procedure to keep the equilibrium angle 



(2.13) constant at the intersections of the surface with the GBs. 

Consider equations ( |2.2| ), ( |3^ ) and the zero level line of </> (conductor's surface) 
passing through the groove root point a in Fig. [l] (a similar analysis could be performed 
for groove root point b). Since the tangential vector to the zero level line at a (as well as 
to other level lines at x = if <f> is kept a signed distance function), is r = (cos 8q, sm $o)> 
then «, (Ip) imply 



n. 



sin o 



Yj^, n y = cosv 



Dividing the first equation in (|3.5|) by the second one gives 

4>x = ~4> y tan 9 . 



1/2- 



(3-5) 



(3-6) 



As pointed out above, we approximate the spatial derivatives (in normal direction) of <fi 
at the boundaries of the computational domain by second order one-sided differences. 
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Therefore, equation ( |3.6| ) written on the left boundary x = takes the form 

-3^- + y 1J --0 2J = _ ^ +1 -^- l tan m _ x 

2Ax 2Ay 

where m is the number of grid points in the vertical direction, and Ax and Ay are grid 
spacings in the horizontal and vertical directions, respectively. Rearranging the terms 
in ( |3,7|) gives the following set of nonhomogeneous linear algebraic equations with a 
tridiagonal matrix for unknowns (f>oj, j = 0, ...m — 1: 

tan O , 3 tan O , Hij - <h,j to o\ 

~ W^" 1 ~ 2A^°' j + W 0OJ+1 = 2Ax ' (3 ' 8) 

The solution to this (and to the similar set at x = L) is easily and accurately found 
at the beginning of every stage of a Runge-Kutta time marching, thus providing the 
field <p which incorporates the correct equilibrium angles at the groove roots. 
The described procedure allows us to have a straight horizontal line 

y{x, 0) = const. (3.9) 

where const. > gives the initial height of the material domain, as initial condition 
for LS simulations. Note that initial condition (|3.9| ) does not match the boundary 
condition ( [2.13 ). This implies that a singularity exists at x = 0,x = L at t = 0. 



This singularity does not present a barrier in solving the system numerically when we 
select an appropriate numerical scheme. Physically, the equilibrium angle is formed 
instantaneously compared with the time needed for the evolution of the surface. We 
are not concerned with the details of this instance. We could use an initial surface 
which is consistent with the boundary conditions (as done in || where Mullins profile 
served as an initial condition, and we followed the evolution of this profile in time). 
However, the choice of initial condition ( |3.9| ) is more physical. 

To close this section, we present the details of the calculation of the normal velocity 
function (gH) : 

1. Calculate curvature induced flux J^ Ka from ([2.3|). It is nonzero even at the first 
time step, since the equilibrium angles ( 2.13[ ) are formed instantly 



2. Calculate the first term in (3.4) by applying the formula 



dJ7 K 



dr 

f - h" ■ ■ < r. — 2K ..,.<• ' ..< >,. — 



V 



vj7 k -t 



(3.10) 



+ KJj K ' 



K XX 4> 2 + 2K X y<f) X (t)y - Kyy4>l 

^ 2 +( j )2 + K ( K v( n x + n V) ~ K A n V ~ n x)) 

3. Solve the electrical problem, find electrically induced surface flux from ( |2.6| ), 
As pointed out above, the discontinuous distributions of electrical quantities are 
smoothed across the surface profile - by a hyperbolic tangent law: 

font i T{ n Tout Tin . i aJL i \ /oi-iN 

r = 1 tanh/3^(x,y) (3.11) 

where (5 is a large constant adjusting parameter, and r is either the electric in- 
tensity E(x, y) (first electrical model), or the specific conductivity k(x, y) (second 
electrical model) 
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4. Given values of the electrical intensity E along the GBs, calculate electrically in- 
duced GB fluxes, J g b, from ( p. 12 ). Applying boundary condition (2.14), calculate 
corrected values of jf along grid lines x = 0, x = L 

5. Calculate the second term in (EO) by applying the formula (|2.2|) where J~ replaces 



K s . 



4 Numerical results and discussion 

Several comments should be made before we present the results of the numerical sim- 
ulations. 

• Due to the large number of material parameters involved we concentrate on the 
influence of the one which was predicted to greatly affect the grooving process, i.e. 
the ratio of the GB to surface diffusivity, = D g b /D s Q . The parameter set we 
choose for the simulations corresponds to copper, Cm, at temperatures about 600 
K. It should be noted that (i) the experimentally measured values of diffusivities 
could vary, according to different sources, by up to 3 orders of magnitude, and 
(ii) D s can be smaller than D g b, due to, for example, surface contamination, thus 
giving > 1. Accordingly, we fix the value of D s and vary D g b in a wide range, 
thus varying the GB flux J gb (|2.12[) . 

• We study the advancing (elongating) grooves characterized by the positive values 
of J g b (matter flows out of the groove cavity and into the GB). In this case 
the electric field intensity vector is directed upwards (see Fig. ||), the positive 
potential being prescribed on the lower boundary of the computational box and 
the negative potential on the upper one. Reversing the direction of the field 
produces receding grooves or "ridges", Fig. B characterized by a negative GB 
flux (matter flows out of the GB into the groove cavity) . The advancing grooves 
are of more practical interest, as explained in the Introduction. 

• Most of our results are obtained with electrical model 2 (see section |2.1| ), based 
on the solution of Laplace's equation for the potential. As expected, this model 
produced more accurate results than the approximate model 1 based on the 
piecewise-constant electrical field. However, model 1 [] proved to be rather useful 
in our simulations, since it produced qualitatively good results in greatly reduced 
(comparable to the model 2) computational times; for very fine grids (80 x 80 
resolution) the speedup achieved is by as much as a factor of 6. 

Fig. |H (a)-(d) shows the space/time evolution of the initially flat surface of the 

conductor for different values of r^. The parameters are as in Figs. [2], [3| and U + = 

-XJ- = 5.0 x 10" 3 V, k in = 10 8 (Hm)" 1 , k out = 10 7 (IWp 1 , 5 gb = S s = 3.5 x 
10 -iu 

m, z* — z* b — 5. The surface profiles are dumped every 5000 time steps, the 
dimensions of the computational box are 0.5 /im x 0.5 fj,m, and the grid has a 60 x 60 
resolution. 

1 Or even more simplified model of the constant field throughout the entire domain C + O (Fig. [j]) used 
in most, if not all studies of electromigration (see H for example) 
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Figure 4: GB grooving by surface/GB diffusion driven by the surface curvature gradients 
and the electromigration. (a): r d = 0.224, (b): r d = 0.336, (c): r d = 0.561, (d): r d = 22.424. 
The surface profiles are dumped every 5000 time steps. The time labels correspond to the 
(physical) time at which the last profile is dumped. 



In case r d is much less than one (Fig. |4| (a)), i.e. the GB flux is relatively small, and 
we observe that the evolution of the surface is similar to the one shown in Fig. |3[ 
It slows down in time, providing to be not very dangerous in the sense of failure of 
the conductor. The evolution proceeds faster as r d increases (Fig. ||| (b)). After the 
GB grooves merge and form a single profile, this profile starts to advance slowly (note 
that the surface curvatures at the groove roots are still positive, at least during the 
time of the observation). Yet larger values of r d (Fig. || (c)) result in changes of the 
morphology of the surface profile in the near-groove-tip regions. The latter means that 
the sign of the surface curvatures at the groove tips changes from being positive to a 
negative one in a realively short time after the evolution starts, indicating the surface 
tendency to form so-called slits. In the case of Fig. || (c) this transition takes place at 
the time step n* , 15000 < n* > 20000; see also Figures |7[ No qualitative change 
in the surface shapes is observed as r d is further increased - up to the limit where the 
numerical method is applicable (note the significant losses of accuracy in Fig. |] (d), 
which corresponds to r d = 22.424). Fig. || (d) differs from Fig. || (c) only in the 
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increased velocity of the surface's advance and in a more rapid transition from positive 
to negative curvature at the groove roots (the decreased n*). The evolution regime 
shown in Fig. [| (c), (d) is known as the ^-regime 0, ^|. 

In Fig. 5| (a) - (d) we plot the distance d traveled by the groove tip as a function of 
time. Fig. | (a) corresponds to the case of Fig. || (d) (r d = 22.424, L = 0.5 /im, U = 
U + = U~ = 5.0 x 10~ 3 V). One sees that the steady-state velocity of the surface 
advance is attained rather rapidly, within approximately 6 min. Also note that the 
groove tip traveled (in only 1 hour) a distance which is a little less than half the grain 
size. In Fig. [| (b), (c) we investigate the influence of the grain size, L. When r d is 
small and the applied voltage is small too (Fig. | (b), r d = 2.424, U = 1.0 x 10~ 4 V), 
the evolution is driven mostly by the surface flux J^ K . Then, smaller sizes of the 
grains result in larger velocities of the groove tip. This is because the curvature of 
the surface increases as the grain size decreases, resulting in the increase of J^ K ■ The 
example of such a transitive grooving regime (from classical Mullins regime, Fig. || (a) 
to vl-regime, Fig. || (c), (d)) is presented in Fig. [| (b). If electrically induced fluxes jf 

(a) (c) 
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Figure 5: Distance d traveled by the groove tips as a function of time, (a): = 22.424, 
refer to Fig. | (d); (b): r d = 2.424, U = 1.0 x 10~ 4 V; (c): r d = 2.424, U = 1.0 x 10~ 3 V; 
(d): r d = 2.424, L = 0.5 /im. 

and J g b are dominant (^-regime, Fig. [5] (c), r d = 2.424, U = 1.0 x 10~ 3 V) then, in 
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contrast with the previous case, larger grain sizes result in larger groove tip velocities, 
as predicted in Q. The dependence of the groove tip velocity on the applied voltage 
is illustrated by Fig. |5] (d) (r^ = 2.424, L = 0.5 /im). GB grooving proceeds faster 
in strong electric fields due to the amplification of the electromigration and associated 
diffusion fluxes Jf and J g b- 

For completeness, in Fig. |6| (a)-(d) and in Fig. [7| (a)-(d) we present plots of the 
surface curvature, the J^ K and jf diffusion fluxes, and the normal velocity function 
F for the case of Fig. || (c). The data in Fig. |6] correspond to time step 1000 (transient 
stage), while the data in Fig. correspond to time step 35000 (steady-state stage). 
Grooves develop faster during the transient stage (compare Fig. ^ (d) and Fig. |?] (d), 
see also Fig. |5|), when the curvature of the surface in the near-groove tips regions is 
positive and both J^ K > and jf > fluxes tend to elongate the grooves, providing 
the flow of matter out of the groove tips. As the steady-state approaches, the curvature 
of the surface in the near-groove tips regions becomes negative (Fig. [7| (a)) and the Jf 
flux is still into the GB but the flux J^ K changes the direction and slows the evolution 
down (Fig. @ (b), (c)). 
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Figure 6: Plots of the surface's curvature (a), diffusion fluxes (b), (c) and the normal velocity 
function (d). = 0.561, refer to Fig. f| (c). The corresponding time step is 1000. 
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Figure 7: Same as Fig. [| but the data correspond to the time step 35000 at Fig. f| (c). 



Fig. |9] shows the ^4-regime of GB grooving obtained with the use of the electrical 
model 1 (see section 2.1; to be compared with the Fig. || (d)). The computations are 
less accurate if this model is employed, resulting in highly asymmetric surface profiles. 
However, the dynamics of surface evolution could be predicted from these simulations 
and we made a heavy use of the electrical model 1 for trial numerical experiments. It 
is worth to note that the run time to obtain Fig. ^ is 2.2 hours on SGI workstation 
with 194 MHz IP25 processor, compared to 8.1 hours for Fig. (d). 



16 



0.43 
0.42 
0.41 




Figure 8: The ridges formed if the matter flows into the groove tips (r<j = 1.121, the electric 
field intensity vector is directed from top to bottom) . The physical parameters are as in Fig. 

1 




0.45 0.5 



Figure 9: GB grooving, = 22.424 (compare to Fig. |] (d)). The electrical model 1 was 
used. 
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5 Conclusions 



The Level Set method was used to model the GB grooving by surface/GB diffusion in 
an idealized polycrystalline interconnect. The diffusion is driven by surface curvature 
gradients and external applied electric field. The results demonstrate the high potential 
of the LS method for the simulation of complex failure phenomena in microelectronic 
interconnects. The plans for future research are: 

• to obtain more physical results with the current version of the code and compare 
them to experimental ones 

• to improve our numerical procedure to make possible the simulation of the prop- 
agation of slits (i?-regime, [2], The latter are characterized by (almost) 
straight vertical walls [] and negative curvatures at the slit tips. The slits are 
formed at high values of r^, and are supposed to propagate in a local steady 
state, leaving the rest of the surface behind. Physically, the surface of the con- 
ductor cannot accomodate a very large GB flux; the groove tips then become 
diffusively detached from the remaining surface. At the moment our numerical 
procedure does not allow to fully trace the evolution of the slits and the surface 
left behind. In our opinion, a locally refined grid is needed to provide high ac- 
curacy in the near-slit-tip regions; however, the adaptation of the LS method to 
such grids may be not straightforward. In our simulations (which make use of an 
uniform grid), the instability steps in shortly. 

• to incorporate mechanical stress in the analysis 

• to speed up the computations by making use of an implicit scheme for the solution 
of the equation ( |3.1[ ). This will allow larger time steps. 

6 Appendix 

Derivation of the 2-dimensional electrostatic equation 

We consider a conducting strip made of a thin metal film, attached to a strip of 
nonzero conductivity substrate. The metal film may be continuous or it may be made 
of conducting patches with voids in between. We allow the metal film and substrate to 
have variable thickness. In the present formulation we neglect the interface resistance. 
The electrodes are attached to the strip and to the substrate. We may want to compute 
the local field strength which determines the resulting electromigration. This is a 
more realistic model then the model based on the assumption of a zero conductivity 
substrate. It allows us also to consider the behaviour of a metal film with varying 
effective thickness at no extra cost. 

6.1 The 3-dimensional problem 

The 3-dimensional Problem Ohm's law implies: j = crE = — aV3<p, where j is the 
electric current density vector, E is the electric field vector, <f> is the electric potential 

2 The flux jJ K is zero along these walls. 
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and a is the material conductivity. For steady fields Maxwell's equations with vanishing 
space charge give: 

V,.?=0,»ta.V, = («) 

Hence 

V 3 • (<rV 3 </>) = 0. (6.2) 

At all external (lateral) boundaries there is no flux in the direction of the normal, n, 
so that n ■ j = 0, and using ( |6.1| ) one gets: 



n-V 3 cp = 0. (6.3) 

The conditions Q6.3| ), together with values of the potential specified at the strip and 
substrate edges and the continuity and jump conditions at the interface, constitute 
boundary conditions for equation ( |6.2D in the two layers. Thus the three dimensional 
potential can be found, in principle, as the solution of a well-posed three-dimensional 
boundary value problem. However such a solution can be very expensive to get in 
the present geometry, in particular as singularities in the solution will appear at sharp 
geometrical corners at crystal boundaries or voids, requiring high resolution or com- 
plicated integration formulae. To avoid this (probably unrealistic) behaviour of the 
solution and to avoid solving three dimensional problems many times, as required by 
the time development of the process, we proceed with an approximate approach sug- 
gested by (singular) perturbation theory. 

6.2 The 2-dimensional equation 

We assume that <f> and a change over a characteristic length scale L in the horizontal 
directions x and y but over a scale H in the vertical. Furthermore we assume that 
e = H/L <C 1. Using scaled variables in (|6.2[), 



we get: 



(X,Y,Z) = (x/L,y/L,z/H) (6.4) 



e V 2 (.V 2 « + A ( ff *) _ o, where V 2 = ( A, JL) . (M) 
Singular perturbation analysis considers an expansion 

(j)= 0o + £ 2 </>i + e 4 </>2 + - (6.6) 



where (fik are functions of order 0(1) in e. Substitution of ( |6.6| ) in ( |6.5| ) gives relations 
for the functions (fik by grouping terms according to their order in e and equating each 
group to zero. The zeroth order term gives : d 2 <fto/dZ 2 = 0, thus 4>q is a linear function 
in z for every x and y while taking into account (|6.3j ) kills off the z dependence, so 
that: 

( j )0 = MX,Y). (6.7) 

Thus at this stage 4>o is an arbitrary function of the horizontal coordinates X and Y . 
The first order equation and the boundary conditions in Z result ultimately in the 
approximate two-dimensional equation for 4>q Jn 



V 2 (/ii<ri + h 2 (J2)V2(po = (6.8) 
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where h\,a\ and h2,o~2 are ) respectively, the heights and conductivities of the two 
layers under consideration. The equation ( |6.£| ) is solved with boundary conditions in 
the (X, Y) plane. We remark that the approximate independence of the potential cj> on 
the Z coordinate justifies also the two dimensional approach for the electromigration 
equation. This behaviour is a consequence of the small aspect ratio assumption and 
the normal derivative boundary conditions Q6,3j ), where one must also involve a small 
slope assumption. 
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